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Abstract 

Numerical computations of stationary states of fast-rotating Bose-Einstein condensates require 
high spatial resolution due to the presence of a large number of quantized vortices. In this pa- 
per we propose a low-order finite element method with mesh adaptivity by metric control, as an 
alternative approach to the commonly used high order (finite difference or spectral) approxima- 
tion methods. The mesh adaptivity is used with two different numerical algorithms to compute 
stationary vortex states: an imaginary time propagation method and a Sobolev gradient descent 
method. We first address the basic issue of the choice of the variable used to compute new metrics 
for the mesh adaptivity and show that refinement using simultaneously the real and imaginary 
part of the solution is successful. Mesh refinement using only the modulus of the solution as 
adaptivity variable fails for complicated test cases. Then we suggest an optimized algorithm for 
adapting the mesh during the evolution of the solution towards the equilibrium state. Consider- 
able computational time saving is obtained compared to uniform mesh computations. The new 
method is applied to compute difficult cases relevant for physical experiments (large nonlinear 
interaction constant and high rotation rates). 
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1. Introduction 

Recent research efforts in the field of condensed matter physics were devoted to the study of 
quantized vortices nucleated in a Bose-Einstein condensate (BEC). Several groups have produced 
vortices in different experimental set-ups flu a is a, leading to numerous theoretical and 
numerical studies aimed at a better understanding of such macroscopic superfluid systems with 
quantized vorticity. 

A typical experimental BEC configuration with quantized vortices is the rotating condensate. 
The condensate is confined by a magnetic potential and set into rotation using a laser beam, 
which can be assimilated to a spoon stirring a cup of tea. Since the solid body rotation is not 
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possible in a superfluid system, the condensate has the choice between staying at rest and rotating 
by nucleating quantized vortices. The number and shape of vortices depend on the rotational 
frequency and the geometry of the trap. The fast rotation regime is particularly interesting to 
explore since a rich variety of scenarios are theoretically predicted: formation of giant (multi- 
quantum) vortices, vortex lattice melting or quantum Hall effects. This regime is experimentally 
delicate to investigate 10,13,0], making numerical simulations very appealing in depicting vortex 
configurations for fast rotations. 

However, numerical simulations of fast rotating condensates are also very challenging for at 
least two reasons. The first difficulty comes from the presence in a condensate of a large num- 
ber of vortices when high rotation frequencies are reached. An example of such configuration 
is illustrated in Fig. Q]for a condensate trapped in a harmonic magnetic potential. We recall 
that a quantized vortex is a topological defect of the macroscopic wave function describing the 
condensate: 

<A = ^p(x,y,z)e m ^\ (1) 

where p is the local atomic density and 9 the phase. In other words, p = in the core of the vortex 
(no condensed atoms are present) and around the vortex there exists a frictionless superfluid flow 
with a discontinuous phase field. As a consequence of this phase discontinuity, the circulation 
around a vortex is quantized 

r= £\.dl = n-, (2) 
J m 

where v = i^V(9 is the local velocity (defined by analogy with classical fluids), h is Planck's 
constant, m the atomic mass and n an integer (the winding number). A numerical system has to 
offer sufficient spatial resolution, not only to capture the large gradients of the density p in the 
small-size vortex cores, but also to cope with phase discontinuities that extend up to the edge 
of the condensate (see Fig. Q]). This explains the use in the literature of discretization methods 
with hi gh o rder spatial accuracy: Fourier spectral 11211 . sixth-order finite differences 



ccuracy: 
ii, 



1 1311141 Hal, sine-spectral IU61I17H . Laguerre-Hermite pseudo-spectral 111811 . etc. 

The second numerical difficulty in computing such configurations comes from the numerical 
algorithm used to converge to stable states with vortices. Most of the numerical algorithms 
proposed in the literature use the so-called imaginary time propagation of the wave function. 
A typical computation (Fig. |2]i starts from an ad-hoc initial configuration and iteratively search 
for a minimizer of the energy describing the system (such methods are described in the next 
section). During the iterative process, the vortices move slowly in the condensate towards their 
final equilibrium locations. Depending on the initial condition, new vortices could also enter 
the condensate. This is the case in Fig. [2] where a converged computation for a lower value of 
the nonlinear interaction constant is used as initial condition. This evolution, called imaginary 
time evolution since it has no physical relevance, has to be accurately captured by the numerical 
system and brings up the question of the behavior of standard dynamic mesh adaptivity methods 
in this context. To the best of our knowledge, this question was not addressed in the literature. 

In this paper we tackle the two above mentioned difficulties by using a low-order finite el- 
ement method with mesh adaptivity, as an alternative of commonly used high-order methods. 



Finite element method have been already used II 19112011 to compute vortex states in rotating BEC, 



but with fixed meshes. An attempt to adapt the mesh was made in 12111 by using a fixed computa- 
tional domain with different mesh densities; finer meshes were initially set in subdomains where 
vortices are guessed to lie in the final equilibrium configuration. 
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Figure 1: Example of fast rotating condensate (harmonic trapping potential, g = 5000, Q./ai ± = 0.95) computed with the 
present method. Contours of atomic density p (left, low density in black) and phase 8 (right) of the converged (stationary) 
state. Note the dense Abrikosov vortex lattice and phase discontinuities joining the border of the condensate. 
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Figure 2: Illustration of the imaginary time evolution of the solution before reaching the converged state displayed in 
Fig. [T] Energy decrease and contours of atomic density p for intermediate states. Note the nucleation of new vortices 
and their rearrangement in a more and more regular Abrikosov lattice. The jumps in the energy evolution correspond to 
mesh refinement. 
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It is important to note that the mesh adaptivity is also of great interest for the simulation of 
vortex states in type-II superconductors. Such systems are described by the Ginzburg-Landau 
(GL) macroscopic model that has close resemblance with the GP equation when high values 
of the GL parameter (kappa) are considered 02211 . Several key studies 123L 24] have set the 
mathematical and algorithmic basis for the use of finite element method to simulate vortex con- 
figurations governed by the GL model (see |25] for a review). However, as mentioned in l25ll . 
using mesh adaptivity in computing vortex lattices with a large number of vortices is still a com- 
putational challenge in this field too. 

The purpose of the present approach is to use a dynamic mesh adaptivity that allows to 
follow the evolution of vortices during the computation until convergence. To this end, we start 
by implementing in a low-order (piecewise linear) finite element setting two different algorithms 
to compute stationary vortex states: a classical method based on the imaginary-time propagation 
of the wave function and a Sobolev gradient descent method for the direct minimization of the 
energy functional. These two algorithms are described in the next section. Section[3]presents the 
finite element setting and the mesh adaptivity strategy based on metric control. Several numerical 
experiments are designed in section [4] We start by answering the basic question of the choice 
of the variable used to adapt the mesh. In particular, we show that the approach, that might 
appear as natural, of refining the mesh following the atomic density p is not always successful. 
Extensive numerical tests prove that refinement using simultaneously the real and imaginary part 
of the solution as adaptivity variables is the successful approach. The new adaptive mesh strategy 
is shown to bring an important computational time saving over computations using refined fixed 
meshes. Finally, the proposed method is applied to compute difficult cases, with large nonlinear 
interaction constant and high rotation rates, that are relevant for physical experiments. 



2. Numerical methods to compute minimizers of the Gross-Pitaevskii energy 

2.1. Mathematical problem 

In the zero-temperature limit, a dilute gaseous BEC is mathematically described by a macro- 
scopic complex wave function iff(x), which spatial configuration is obtained by minimizing the 
Gross-Pitaevskii (GP) energy. We consider a BEC of N atoms trapped in a magnetic potential 
Vtmp with radial symmetry and transverse trapping frequency ai ± . The condensate is rotating 
along the z-axis with the angular velocity Q.. It is common practice to scale the variables using 

as characteristic length the harmonic-oscillator length d = J— —, where h is Planck's constant 

° Y mcox 

and m the atomic mass of the gas. Using the scaling, r = x/d, u(r) - tf/(x)d 3/2 / yN, Q = Cl/a> ± , 
we obtain the non-dimensional energy (per particle) in the rotating frame: 

E(u)= { )-\Vu\ 2 + V lrap \u\ 2 + ^\u\ A -£liu*{A<V)u, (3) 
Jo 1 1 

where V tmp = j^-V, rap , and A - (y, -x, 0). We denote by u* the complex conjugate of u. The 
interactions between atoms are described by the constant g = 47r ^"' , with a s the s-wave scattering 
length. The mass conservation constraint becomes in this scaling: 



r m 2 =m 2 =i, (4) 

Jo 
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where we denote by ||.| = ||.||l 2 (:d,C)- Note that we have considered that u(r) — > 0, as r — > oo and, 
consequently, the condensate could be confined in a bounded domain T>. 

We consider in the following the two-dimensional problem defined on T) c K 2 , with ho- 
mogeneous Dirichlet boundary conditions u = on dD. For given constants Q., g and trapping 
potential function V trap , the minimizer u g of the functional (O under the constraint is called 
the ground state of the condensate. Local minima of the energy functional with energies larger 
that E(u g ) are called excited (or metastable) states of the condensate. 

We present in the following two different methods to compute minimizers of the GP energy. 

2.2. Imaginary time propagation: Runge-Kutta-Crank-Nicolson scheme 

Most of the numerical algorithms proposed in the literature to compute minimizers of the 
GP energy use the so-called normalized gradient flow lfl6ll . It consists in applying the steepest 
descent method for the unconstrained problem, 

du 1 dE(u) V 2 w 9 . 

— = ---T^ = — - VtrapU - g\u\ 2 u + iQA'Vm, (5) 
at 2 ou 2 

to advance the solution u e C from the discrete time level t„ to f„+i; the obtained predictor 
u(r, f„ + i) is then normalized and used to set the solution at t n+ \ satisfying the unitary norm con- 
straint: 

, . a u(r, t„ + \) 
u(r,t n+ {) = — — — -. (6) 
IKr,? n+ i)|| 

It is interesting to note that <(5j is commonly referred as the imaginary time evolution equation, 
since the right-hand side corresponds to the stationary Gross-Pitaevskii equation. The gradient 
flow equation (O (or the related continuous gradient flow equation, see lfl6ll ) can be viewed 
as a heat equation in complex variables and, consequently, solved by different classical time 
integration schemes (Runge-Kutta-Fehlberg iflolrTlll . backward Euler lfl6l JjJ 17, 18], second- 



order Strang time-splitting lfl6l [l9ll . etc.). We describe in the follo wing the combined Runge 



Kutta-Crank-Nicolson scheme that was successfully used in II13L 1141 115H to compute stationary 
three-dimensional BEC configurations for different trapping potentials. 
If we write (f5]l under the general form 

— =N(u) + £(u), (7) 
ot 

with N{u) containing non-linear terms and _£(«) linear terms, a combined three-step Runge-Kutta 
and Crank-Nicolson scheme reads: i26[ l27ll : 

Uk+\ — Ufc Cfc 

= a k N(u k ) + b k N(u k -\) + —-C(u k+ i + u k ), (8) 

ot v ' 2 

Runge-Kutta 



Crank-Nicolson 



where k — 1, 2, 3 are the substeps needed to advance the solution from t„ to t n+ \. The following 
values for the coefficients 



ensure the third-order accuracy in time for the Runge-Kutta part and second-order overall accu- 
racy. Note that the intermediate integration time values are t k = t n + c k 6t, with c\ + c% + C3 = 1. 
An important computational advantage is that the scheme is low-storage and self-starting. In- 
deed, since b\ — the storage of the solution u"~ l at the previous time-step is not necessary. For 
numerical purposes, the equation to solve is written as: 

1 Q 



2 -Cj q k = [a k N(u k ) + b k N(u k -\)~\ + c k £(u n ), q k = u k+l - u k , (12) 

with the variational formulation: find q k e H l Q {D, C) such that Vv e H l Q (D, C), 

f ^qkV - ^-L{qk)v = f (a k N(u k ) + b k N(u k -i))v + c k £(u k )v. (13) 
Jv[ot I J Jo 

Depending on the choice of the linear operator in (0, we can distinguish between different 
schemes. In lfl3l IT4L [l5ll the linear operator was defined in the classical way: £(u) — V 2 (m). 
We use in the following a different choice that resulted in a better stability of the scheme: 

£(«) = V 2 (m) + 2/QA'Vh, (14) 

N(U) = -2 [g\u\ 2 + V trap ] U. (15) 

For this method, the mass conservation constraint © is taken into account by using the discrete 
normalization ©. 

2.3. Direct minimization: Sobolev gradient descent method 

Another method to compute stationary BEC states is to directly minimize the GP energy 
(0 using steepest descent methods. It is interesting to note that in the descent method (0), the 
right-hand side represents the L 2 -gradient (or ordinary gradient) of the energy functional. An 
important improvement of the convergence rate of the descent method is obtained by replacing 
the ordinary gradient with the gradient defined on the Sobolev space H l (D, C). The reason 
is that the use of Sobolev gradients is equivalent to a preconditioning of the ordinary gradient 
method. The idea of introducing the Sobolev gradient in a descent method was developed by J. W. 



Neubergerin the 1970's and is now used in several fields of numerical analysis (see B28H ). On the 
related topic of finding minima of the Ginzburg-Landau energy functional for superconductors 
i23ll24ll . the Sobolev gradient method was first presented in ll29tl . Recent developments of the 
method in a finite element setting include the minimization of Schrodinger [ 30] or Ginzburg- 
Landau type functionals Pill . 

In the framework of computing critical points of the Gross-Pitaevskii energy with rotation, 
a descent method based on the H 1 Sobolev gradient was used in lflo[[Tlll . in conjunction with a 
spectral method for the spatial discretization. In [32] we have equipped the Sobolev space H l 
with a new inner scalar-product and used the associated gradient to improve the convergence of 
the descent method for high rotation frequencies. The new inner product is 

(u,v)h a = J (u,v) + (V A u,V A v), (16) 
Jd 

where = V + iQA', and (■, ■} denotes the complex inner product. The new Hilbert space is 
denoted by Ha(D, C). Hence, the gradient of the energy functional satisfies the equation: 

< V Ha E, v > Ha =< V L 2£, v > L 2 . (17) 



The numerical method introduced in 113211 consists in the following steps: 
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We first compute the gradient V ' h a E. Observing that an equivalent definition of the Ha 
scalar product is: 



u,v> Ha = \ <[l + £l 2 (y 2 + x 2 )]u,v) + (Vu,Vv)-2iQ.(A'Vu,v), 
Jo 



(18) 



we infer that the gradient @ = Vh a E could be directly computed from (JT7J as the solution 
of the variational problem: 

f [l + Q 2 (y 2 + jr)l Qv + VQVv - 2i£l(A'Vg)v = RHS, Vv e H^D, C), (19) 
jd 

where the right-hand-side term represents the L 2 gradient (in the weak form): 
RHS = I VwVv + 2 \Vtrap u + (g\u\ 2 )u - iCM'Vu] v. 

JD 



(20) 



• In order to satisfy the mass constraint (0), we project the gradient Vh a E onto the tangent 
space associated to the constraint. In our case, we project onto the null space of /?'((<)> 
where (5{u) = J D \u\ 2 . The final expression (see 13211 for details) of the projection that will 
be used for numerical implementation is: 

P u ,h a = G- ^ v Ha , (21) 

&(u, v Ha ) L 2 

with % denoting the real part, and Vh a computed such as that 

K(v Ha , v) Ha = p'{u)v = %{u, v) L 2 . (22) 

• The solution is finally advanced following the general descent method: 

=u n -6tP u , HA 0(u n ). (23) 

It should be noted that the projection step ensures that the norm of the initial condition (u°) is 
preserved through the iterative process ( l23l . 



3. Finite element spatial discretization and mesh adaptivity 

The finite element implementation uses the free software FreeFem++ l33ll . which proposes 
a large variety of triangular finite elements (linear and quadratic Lagrangian elements, discon- 
tinuous P 1 , Raviart-Thomas elements, etc.) to solve partial differential equations (PDE) in two 
dimensions (2D). FreeFem++ is an integrated product with its own high level programming lan- 
guage with a syntax close to mathematical formulations. FreeFem++ was recently used to test 



algorithms for the minimization of Schrodinger or Ginzburg-Landau functionals 130LI31I1 



7 



3.1. FreeFem++ implementation 

It is very easy to implement the variational formulations associated to the above described 
algorithms using FreeFem++. We outline here the main features of the finite element implemen- 
tation that were helpful in writing efficient FreeFem++ scripts. Let T/, be a family of triangula- 



tions of the domain D. We assume that 7), is a regular family in the sense of Ciarlet [34], with 
h > belonging to a generalized sequence converging to zero. We denote by P'(T) the space of 
polynomial functions on triangles T e T/,, of degree not exceeding I > 1 . We also introduce the 
finite element approximation spaces: 

W[ = {w h e C°(D h ); w h \ T € P'(T), Vr e T h ) , (24) 

and 

V l h = {w h €W l h ;w h \r h =o). (25) 

The finite dimensional space V' h is a subspace of H^{D) and therefore will be used to discretize 
the variational formulations previously written. We use in the following P 1 (I = 1, piecewise 
linear) finite elements to approximate the solution and a P A representation of the nonlinear terms. 
It is interesting to note that FreeFem++ allows to switch to P 2 (Z = 2, piecewise quadratic) finite 
elements by a simple change of the definition of the generic finite-element space W' h . 

An efficient implementation of the algorithms described in the previous section is obtained 
using the pre-computation of the complex matrices associated to linear systems. For the imagi- 
nary time propagation method, the integral form ( fT3l leads to the following linear system: 



1 A C * A C k A 

St + ~2 ~ 2~ 



Q k = A 4 M .(a k N k + b k N k -i) - c k A c U k + c k A n U k , (26) 



where U k is the solution vector at substep k of the Runge-Kutta method and Q k = U k+ \ - U k . 
Denoting by w\ the basis functions of the space Vl, the matrices in d26l i are defined in the classical 
way using 1=1: 



lc)m,p = V(w>) m V(w],V (28) 



(Am),,,,, = f (w l h ) m (wl) p , (27) 

•JO,, 

(An) m , P = (2/Q) f {A'V{ w l h ) p ){ w \) m . (29) 

Jo,, 

Nonlinear terms N k , corresponding to (1151 1. are computed with higher accuracy using P 4 finite 
elements. The (non squared) matrix A 4 M is consequently computed as: 

(A 4 M ) m , p = f (wl) m (w 4 h ) p . (30) 
Jo, 

Previous two-dimensional integrals are computed using a fifth order quadrature formula. If the 
imaginary time advancement is conducted with fixed size time step, a further optimization comes 
from the storage of the three matrices of the linear systems corresponding to each substep of the 
Runge-Kutta integration procedure. 
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For the Sobolev gradient method, the discrete form of $19[ becomes: 



A S G= A A M .N n + A C U„ - A n U„, (31) 

with N„ corresponding to a P 4 representation of nonlinear terms 2 (v tra p + g|M„| 2 ) u„. The matrix 
As of the linear system: 

(As W = f [l + Q V + x 2 )l (w{) m (w\) p + V(w> ) m V(w> )„ - 2/Q(A'V)(w>) p ( w>) m , (32) 

is computed by a fifth order quadrature formula. An important computational tine saving is 
obtained if the matrix As is stored and factorized before the time loop ( f23l . 

The last point to emphasize concerning the FreeFem++ implementation is that all previous 
equations are solved in complex variables. As a consequence, the corresponding matrices also 
have complex elements. The approach used in J30ll . based on the separation of the real and imag- 
inary part of the unknown variable, results in considerably larger computational times. Besides, 
this separation is not possible when computing the Ha gradient from ( TT9l i. 



3.2. Adaptive mesh refinement strategy 

Mesh adaptivity by metric control is a standard function offered by FreeFem++. Details on 
the ingredients used in the metric mesh adaptation can 

be found in JHSIH HUSH- The 



key idea is to modify the scalar product used in an automatic mesh generator to evaluate distance 
and volume, in order to construct equilateral elements according to a new adequate metric. The 
scalar product is based on the evaluation of the Hessian H of the variables of the problem. 
Indeed, for a P l discretization of a variable^-, the interpolation error is bounded by: 

£ = \x - H#lo < c sup sup fH(x)\(y - z).(y - z) (33) 

Terr,, x,y.zeT 

where Yl^x is the P l interpolate of x, is the Hessian of x at point x after being made positive 

definite, and . denotes the dot product. We can infer that, if we generate, using a Delaunay 
procedure (e.g. 03811 ). a mesh with edges close to the unit length in the metric M = j^, the 
interpolation error S will be equally distributed over the edges a, of the mesh. More precisely, 
we have 

\ajMai < 1. (34) 
co 

The previous approach could be generalized for a vector variable ^ = lx\,X2]- After computing 
the metrics All and AI2 for each variable, we define an metric intersection M — M\ D M2, such 
that the unit ball of M is included in the intersection of the two unit balls of metrics M2 and 
Ali. For this purpose, we use the procedure defined in ll39ll . Let A? and v 7 , — 1,2) be the 
eigenvalues and eigenvectors of Aij, j = 1,2. The intersection metric (At) is defined by 

M = . (35) 

where A(i (resp. AI2) has the same eigenvectors than All, (resp. AI2 ) but with eigenvalues 
defined by: 

A\ = maxU, 1 , v? 7 M 2 v\), i=\,2. (36) 
9 



FreeFem++ uses mesh generation tools developed in Q8LI39H with the novelty that the De- 
launay mesh generation procedure introduces an extra criterion to keep the new mesh nodes and 
connectivity maps unchanged as much as possible when the prescribed mesh by the new metric 
is similar to the previous mesh. This reduces the perturbations introduced when the solution is 
embedded by interpolation from the old mesh to the new one. 

The mesh adaptivity strategy used in this work is based on the fact that the energy of the 
solution decreases during the computation to attain a plateau corresponding to a local minima 



(see Fig. [2j». Since we generally use a convergence criterion [13, 14 based on the relative 

change of the energy of the solution, 5E n = (E„ + \ - E„)/E„ < s c , we monitor the same quantity 
to trigger the mesh adaptivity procedure following the next algorithm: 

1. choose a sequence of decreasing values s' > e c , that represent threshold values for the 
mesh adaptivity; 

2. set i = 1; 

3. if e' +I < 6E„ < s' and 6E„ > e c , call the mesh adaptivity procedure; the solution u is 
interpolated on the new mesh and normalized to satisfy the unitary norm constraint; 

4. if step 3 was performed Nad ^ 1 times, increase i to i + 1 . Limiting the number of mesh 
refinements for the same threshold, is necessary since, at step 2, the interpolation on the 
new refined mesh and the normalization of the solution could lead to an increase of the 
value of 5E„. 

As an example, for the computation displayed in Figs. [T]and|2] we fixed the convergence thresh- 
old to e c = 1(T 8 and mesh refinement threshold values to e e {1(T 6 ,5 ■ 1CT 7 ,2.5 ■ 1CT 7 , 1(T 8 }. 
The number of calls for the mesh refinement procedure was N ac / = 3 for each fixed threshold. 
We can notice in Fig. [2]the jump in the energy evolution when the mesh refinement was applied, 
resulting in a faster convergence to the final value of the energy. 

An essential question that remains when defining the mesh refinement procedure is the choice 
of the mesh adaptivity variable x- Since vortices are characterized by small cores in which the 
atomic density rapidly decreases to zero in the vortex center, it may appear obvious to use as 
mesh refinement variable^ = \u\, the modulus of the wave function. We prove by extensive 
numerical tests described in the next section that this approach is not always successful. In 
exchange, the adaptivity strategy considering simultaneously the real and imaginary part of the 
solution to compute the metrics, i.e. x = t M r, «;], proved effective in capturing the right solution 
with an important reduction of the computational time compared to fixed mesh calculations. This 
strategy was applied in computing the complex vortex configuration displayed in Fig. Q] 



4. Numerical experiments 

In computing stationary states of rotating Bose-Einstein condensates, the initial condition mq 



plays a crucial role. It was theoretically proved in 14111 that in a real-time evolution of the rotating 
condensate, the number of vortices attained by the condensate depend upon the rotation history 
of the trap and on the number of vortices present in the condensate initially. This observation 
also holds for the imaginary-time evolution: for the same rotation frequency, different stationary 
states, characterized by closed values of the energy, could be obtained starting from different 
initial conditions. 

Three types of initial conditions are generally used for computing stationary states in a rotat- 
ing BEC: (i) condensate without vortices, with a wave function distribution derived from a phys- 
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ical model, called the Thomas-Fermi approximation; f ii) condensate described by the Thomas- 
Fermi model on which vortices could be artificially superimposed using an mathematical ansatz; 
( Hi) initial state set equal to a converged state for a different rotation frequency Q or a different 
interaction constant g. The computation depicted in Figs. Q]and[2]was performed for g = 5000 
and started from a converged state obtained for g = 2000. 

The Thomas-Fermi approximation consists in neglecting the contribution of the kinetic en- 
ergy in the strong interaction regime (large values of g). A simplified energy functional is ob- 
tained: 

£ TF (P)= f V trap \u\ 2 + f \u\\ (37) 
Jd 1 

with a minimizer corresponding to the so-called Thomas-Fermi atomic density: 

p TO (r) = \uf = (^J^-) . 08) 

where fi is the chemical potential. Since // is a Lagrange multiplier, a relation that allows to 
compute fi is obtained by imposing the mass constraint in ( 1381 ). The initial condition is finally set 
to u (x,y) = yjp TF (x,y). 

This model is also useful in estimating the necessary size of the computational domain. When 
a rotation Q is applied, the Thomas-Fermi approximation d38l ) stands with V trap replaced by: 

QV 

V ef f(r) = V trap (r)-—. (39) 

The resulting radius 7?", corresponding to the point where = 0, is used to estimate the size ro 
of the domain D in simulations {ro > R^) . 

Initial conditions with vortices are obtained b y superimp osing to the Thomas-Fermi wave 



function distribution a simple ansatz for the vortex II13LI14LI15I1 . For example, an initial condition 



with N v vortices of radius e v and centers (x[„y' v ), i= 1, . . . ,7Y V is obtained by imposing 

uo(x,y) = ^p 1F (x,y) ]~~[ u[,(x,y), (40) 



u' v (x,y) = \\0.5 < 1 + tank 



4 

- (n - e v ) 



exp(tf ; ), (41) 



where (r;, 8/) are polar coordinates in the framework centered at {x\ n y' v ). Note that the ansatz is 
written for singly quantized vortices (winding number equal to 1). 

We present in the following different types of numerical experiments. We start with test 
cases reflecting two different imaginary time evolutions: (i) the number of vortices at conver- 
gence remains the same as in the initial condition; (ii) new vortices enter the condensate. These 
experiments will serve to test different strategies for mesh adaptivity and to ascertain the com- 
puting time gain offered by the present method. Finally, the method is used to compute complex 
configurations relevant for physical rotating condensates. 

We also mention that the converged final state is characterized by its energy E(u) and angular 
momentum L : (u) which gives a measure of the rotation: 

L,( M )= f K(iu(A'V)u). (42) 
Jd 
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4.1. Numerical experiment 1 

In laboratory experiments, the condensate is typically confined by a harmonic trapping poten- 
tial Vtrap = r 2 /2. It is easy to see from d39l l that this potential sets a upper bound for the rotation 
frequency, since for Q = 1 the centrifugal force balances the trapping force and the confinement 
of the condensate vanishes. To overcome this limitation, different forms of the trapping potential 
are currently studied, experimentally and theoretically. We use in this experiment a combined 
harmonic -plus-quartic potential lf42lTl~4l [l~5 . 43 1 that allows high rotation frequencies. 

We set the following parameters of the simulation 



g = 500, V, 



trap 



r 2 /2 + r 4 /4, 



Q = 2. 



(43) 



The computational domain is circular of radius R max = 1 .25 ■ , where the Thomas-Fermi radius 
is for this case = 3.4. The initial mesh is generated using M = 200 equally distributed points 
on the border of the domain. 



adapt |U| M=200 
adapt [Ur, Ui] M=200 
no-adapt M=200 
no-adapt M=400 
6th order FD 




150 200 

iterations 



Figure 3: Computation for g = 500, fl = 2 and combined harmonic-plus-quartic trapping potential. Initial condition with 
6 vortices artificially placed at 0.5 ■ R FF . Energy evolution for constant mesh and different adaptive mesh computations; 
the result obtained with a 6th order finite difference method is also plotted for comparison. Density contours (|w|) for 
initial and converged solution (low density in black). 

The computation is depicted in Fig. [3] The initial condition contains an array of 6 singly 
quantized vortices equally distributed on the circle of radius 0.5-7?" . The converged state contains 
the same number of vortices, but with larger cores than initially set, and placed closer to the center 
of the condensate, at 0.33 ■ R^ F . Two computations with fixed mesh (M = 200 and M = 400) 
are run and compared to adaptive mesh computations using as adaptivity variable^- = \u\ and 
X = [u r ,Ui], respectively. Convergence test is set to 6E„ < 2 ■ 10~ 6 for all computations and 
threshold values for mesh refinement are chosen as e e {10~ 2 , 10~ 3 , 10~ 4 , 10~ 5 , 10~ 6 ). Three 
mesh refinements are done for each threshold {Nad - 3). 
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It is interesting to note from Fig. [3] the monotone decrease of the energy which is typical 
for the steepest descent method. This evolution is not affected by the projection method for 
the unitary norm constraint, as showed for the computations using fixed meshes. The mesh 
refinement results in a jump in the energy evolution curve at adaptivity thresholds. As already 
stated, this is the consequence of the interpolation on the new refined mesh and the normalization 
of the interpolated solution. Such jumps are naturally less visible close to the convergence, when 
small variations of the energy are monitored. 

In order to assess for the correct behavior of the numerical system, we also compare present 
finite element results with those obtained using a high order finite difference method. For this 
purpose, the imaginary time -propagation method presented in section lZ2l was implemented using 
for the spatial discretization a 6th order compact finite difference scheme that offers spectral-like 
accuracy [44]. The method has similarities with that used in lfl3l Ffl [l5ll to compute stationary 
vortex states in a three-dimensional BEC. The finite difference method uses a squared computa- 
tional domain of size 2R max and a uniform mesh of 105 x 105 grid points. The constant mesh 
size Sx - 6y - 0.08 thus becomes similar to the minimum edge size of the final refined finite 
element grid (h mm = 0.08). 

All computations lead to identical configurations of the final, converged state, as represented 
in Fig. |3] A detailed comparison between finite element and finite difference results is offered in 
Fig. [4] The finite element grid contains initially 7054 triangles and ends with an adapted mesh 
with 3722 triangles, while the finite difference mesh has a fixed size of 11025 grid points. A 
zoom inside the zone containing two neighboring vortices of the final configuration shows that 
contours of the atomic density \u\ are almost identical. It should be noted that in adapting the 
finite element mesh, one could impose the values for h„ mx and h mm , which are the maximum and, 
respectively, the minimum edge size of the triangular mesh. Reducing the value of h max will 
result in a finer mesh and smoother contour lines, comparable to those obtained with the high 
order finite difference discretization. However, the present comparison is more than satisfactory 
with a final finite element grid using almost 3 times less grid points than the finite difference 
setting. 




Figure 4: Computational case depicted in Fig. fj] Comparison between the results obtained with the finite element 
method for M = 200 and mesh adaptivity using x — l u r, u{\ and the 6th order finite difference method with a 105 X 105 
equally spaced grid. Details of the contours of the atomic density \u\ and corresponding grids (dashed lines for the finite 
difference results). 
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Run case 




Sobolev gradient method 


Imaj 


binary time method 


M 


N, 


E(u) 


L z {u) 


iter 


CPU 


E(u) 


L z (u) 


iter 


CPU 


adapt [u r , Ui] 


200 


3722 


11.87 


5.118 


232 


55 


11.87 


5.112 


139 


54 


adapt [\u\] 


200 


2586 


12.04 


5.095 


241 


44 


12.02 


5.088 


142 


40 


no-adapt 


200 


7054 


11.98 


5.126 


223 


72 


11.91 


5.085 


75 


43 


no- adapt 


400 


27674 


11.91 


5.169 


243 


315 


11.83 


5.125 


92 


211 



Table 1: fl = 2: run cases corresponding to the numerical experiment depicted in Fig. [3] Parameters of the initial 
mesh (number of points M placed on the border of the circular domain to generate the mesh and number of triangles 
Nt), energy E(u) and angular momentum L z (u) of the final state, and computational efficiency (number of iterations and 
computational CPU time). 



The exact values of the energy E(u) and angular momentum L z {u) characterizing the final 
state are shown in Tab. Q] Compared to the fixed mesh computation using a refined mesh 
(M = 400), the adaptivity strategy using two variables (x = [u r , u[\) gives the closest energy 
value. We can also see from Fig. [3]that this is also the case when comparing with the 6th order 
finite difference result. Meanwhile, this adaptive mesh strategy results in a computational time 
reduction by a factor of 6 for the Sobolev gradient method and by a factor of 4 for the imaginary 
time propagation method. Table Q] also shows that the two numerical methods used to compute 
stationary states behave similarly. Since this is also the case for all subsequent numerical exper- 
iments, we discuss in the following, for the sake of simplicity, only the results obtained with the 
Sobolev gradient method. This method has also the advantage to allow a constant time step for 
different mesh densities (see also II32I1 ). 

The mesh evolution for the two adaptivity strategies can be followed in Fig. [5] Only meshes 
for the first (e = 10~ 2 ) and final (e = 10 ) thresholds are represented. It can be easily seen 
that adaptivity taking into account only the modulus of the solution results in a very dense mesh 
in the center of vortices. Adaptivity following the real and imaginary part of the solution also 
generates a refined mesh in the core of vortices, but also a dense mesh from vortices towards the 
border of the condensate. This allows to have a better representation of the phase of the solution 
(as previously pointed out when discussing Fig. [TJ. We shall see in the following that this feature 
is crucial for the success of the adaptivity strategy when more complicated cases are computed. 
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Figure 5: Mesh evolution during the computation for experiment 1 (see Fig. [5)- First (e = 10 2 ) and final (e = 10 5 ) 
refined meshes are represented for the adaptive mesh strategy using^ = \u\ (a and b) and^ = [u r , «;] (c and d). 

4.2. Numerical experiment 2 

In this experiment, we consider the same parameters as for experiment 1 and increase the 
rotation frequency to Q = 2.5. The initial condition is the converged state previously computed 
for SI = 2. For this case, new vortices are nucleated inside the condensate and the final state 
contains a second circle of 10 vortices. Figure [6] shows that only the adaptive mesh strategy 
based on x — [u r , Ui] converges to a similar stationary state as the computation using the fixed 
refined mesh (M = 400). This is also visible from Tab. [2] when comparing the values of the 
energy and angular momentum of the final state. In exchange, mesh refinement using x — M 
do not allow the nucleation of new vortices; as a consequence, the energy of the system is not 
decreasing and the final state has the same configuration as the initial condition. It is important 
to note from Tab. [2]that the successful adaptive mesh strategy allows a tremendous (factor of 10) 
gain of computational time. 

The explanation for the failure of the adaptive method based solely on the modulus of the 
solution is offered in Fig. [7J A computation is subject to inherent numerical perturbations that 
will trigger the nucleation of new vortices. Such perturbations usually have small amplitudes, 
and the refinement based on the modulus of the solution will damp them since the mesh size in 
these regions is not small enough to capture them. The adaptive mesh strategy using x - [u r , w,] 
generates refined meshes over larger regions than the core of vortices (see Figs. [7}; and[7ji) and 
consequently allow the nucleation of new vortices. 
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iterations 

Figure 6: Computation for g = 500, fl = 2.5 and combined harmonic-plus-quartic trapping potential. Computations start 
from the converged state obtained for fl = 2. Energy evolution for constant mesh (M = 400, dashed line) and different 
adaptive mesh computations: x — M (dash-dot line) and^ = [u r ,u{\ (solid line). Density contours (|«|) for initial and 
converged solution (low density in black). 

An intriguing question that one could raise after analyzing numerical experiments 1 and 2 
is whether the adaptive mesh strategy based on the modulus is successful if the perturbation 
necessary to nucleate vortices are present in the initial condition. This question is addressed by 
performing computations starting from an initial condition with three arrays containing 6, 12, 
and 36 vortices, respectively. The external circle of vortices plays the role of a dense perturba- 
tion field that could trigger vortices for this high rotation frequency. Figure [8] shows that, once 
again, only the adaptive mesh strategy considering simultaneously the real and imaginary pert 
of the solution is successful. The converged configuration for this computation is very similar 
to that obtained when using a refined (M = 400, h mm = 0.0506) fixed mesh or a 6th order finite 
difference method using a 125 x 125 uniformly spaced grid (Sx = 8y = 0.053). 



Run case 




Initial condition 1 


Initial condition 2 


M 


Nt 


E(u) 


L z (u) 


iter 


CPU 


E(u) 


L z (u) 


iter 


CPU 


adapt [u r , 


200 


8968 


6.08 


11.95 


1266 


321 


5.94 


12.87 


222 


195 


adapt [\u\] 


200 


2540 


9.43 


5.41 


456 


33 


7.14 


11.30 


3280 


947 


no-adapt 


400 


27654 


6.23 


12.81 


3041 


3368 


6.31 


13.01 


249 


327 



Table 2: fl = 2.5. Same legend as for tablefT] Initial condition 1 is the converged state obtained for fl = 2 (Fig. |6) and 
initial condition 2 contains an artificially generated state with 3 arrays of vortices (Fig. [8}. 
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Figure 7: Q. = 2.5 Mesh evolution during the computation for experiment 2 (see Fig. [6}. Refined meshes corresponding 
to thresholds e = 10~ 3 and e = 10 . Adaptive mesh strategy using^ = \u\ (a and b) and^ = [u r , «,] (c and d). 



4.3. Condensates with giant vortex or dense vortex lattice 

In order to assess for the efficiency of our numerical system, we consider in this section 
two cases closer to experimental configurations. Such cases are difficult to compute since they 
involve high values for the atomic interaction constant g and/or rotation frequency Q. 

The first case considers the condensate trapped in the harmonic -plus-quartic potential d43l >. 
but with higher atomic interaction constant, g = 1000. Figure [9] shows the evolution of the sta- 
tionary state of the condensate when the rotation frequency is increased. Vortices in the center 
of the condensate progressively merge to form a giant hole, also called giant vortex. This in- 
triguing configuration has been intensively studied in the physical literature 14211 14 1 15l 14311 . The 
adaptive mesh refinement is very useful in computing such cases since the atomic density in a 
large zone in the center of the condensate is close to zero. As a consequence, large triangles are 
generated in the center of the condensate, while the mesh is highly refined in the annulus zone, 
where vortices nucleate. For instance, the simulation for Q = 4 started with an initial mesh with 
N, — 18 670 triangles and ended with a fine mesh with N, = 69 859 triangles. For reference, a 
constant mesh that offers a similar mesh density in the annular zone is obtained for M = 600 and 
contains N, = 108 212 triangles, since all the computational domain is finely meshed. 

The second configuration considers the case, displayed in Figs. Q]and[2] of the condensate 
trapped in harmonic potential and rotating at Q. = 0.95. We recall that for this case the rotation 
frequency cannot exceed O = 1 . The difficulty for this case is to increase the atomic interac- 
tion constant g that sets the amplitude of the nonlinear term. Figure [10] displays the converged 
configurations for increasing g = 5000, 10000 and 15000. The condensate becomes larger with 
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increasing g, and, consequently, contains more and more vortices that arrange into a regular tri- 
angular lattice (or Abrikosov lattice). The large number of vortices present in the condensate 
requires refined meshes making the computations very costly. For reference, the final refined 
meshes contain, for the three cases, N, = 238 262, 405 405, and, 620706 triangles, respectively. 
Nevertheless, such computations performed with FreeFem++ remain affordable on a single pro- 
cessor computer. 

5. Summary 

We have shown in this work that low-order finite element methods with mesh adaptivity 
are a valid alternative of commonly used high-order methods in computing stationary vortex 
states of a fast-rotating Bose-Einstein condensate. The mesh refinement using metric control 
proved effective in computing difficult cases with a large number of vortices or with giant vortex. 
We showed by extensive numerical tests that adaptive mesh strategy using simultaneously the 
real and imaginary part of the solution to compute metrics is the successful approach. The 
strategy based only on the modulus of the solution failed for complicated test cases. An effective 
algorithm for mesh adaptivity was proposed, with an important computational time reduction 
over computations using refined fixed meshes. 

The present finite element discretization with mesh adaptivity was tested with two numerical 
methods for computing stationary states: a Sobolev gradient descent method for direct minimiza- 
tion of the energy functional and a method based on the imaginary time propagation of the wave 
function describing the condensate. The method is, however, of more general interest, and could 
be used in conjunction with different numerical methods for computing imaginary or real time 
evolution of superfluid systems with vortices, such as rotating Bose-Einstein condensates or type 
II superconductors. In this context, it is interesting to mention that, after the present manuscript 
had been completed, the recent review paper lE5ll was brought to our attention. Among the re- 
maining issues in developing numerical methods for computing vortex states in superconductors, 
adaptive methods were considered in 02511 as challenging because of the complicated patterns of 
the solution with vortices. The necessity to refine the mesh not only around vortex cores was 
intuitively recalled when discussing the different patterns displayed by the the real and imagi- 
nary parts of the solution. The present study confirms in some sense this intuition and offers 
an effective method to answer the challenging question of computing solutions with quantized 
vortices. 
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Figure 8: Computations with the same parameters as in Fig. [6] but starting now from an artificial initial condition 
containing three arrays of vortices. Only the adaptivity strategy with^ = [u r ,uf\ allows to obtain a correct final state, 
which is almost identical to the configuration obtained with a refined fixed mesh (M = 400) or a 6th order finite difference 
(FD) method with a 125 X 125 uniform grid. 
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Figure 9: Condensate trapped in a harmonic-plus-quartic potential (g = 1000). Two- and three-dimensional represen- 
tation of the atomic density contours (low density in black) for increasing values of the rotation frequency CI. Note the 
formation of a giant vortex (hole) in the center of the condensate. 




Figure 10: Condensate trapped in a harmonic potential (Q. = 0.95). Two- and three-dimensional representation of the 
atomic density contours (low density in black) for increasing values of the atomic interaction constant g. Note the 
increase of the number of vortices with increasing values of g. 
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